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A Monte Carlo method is used in order to simulate the 
competition between the molecular relaxation and crystal- 
lization times in the formation of a glass. The results show 
that nucleation is avoided during supercooling and produce 
self-organization in the sense of the rigidity theory, where 
the number of geometrical constraints due to bonding and 
excluded volume are compared with the degress of freedom 
available to the system. Following this idea, glass transitions 
were obtained by producing self-organization, and in the case 
of geometrical frustration, self-organization is naturally ob- 
served. 



I. INTRODUCTION 

When a liquid melt is cooled, usually two things can hap- 
pen: the melt crystallize, or if the speed of cooling is high 
enough to avoid equilibrium, a solid without long range or- 
der (a glass) is formed. This last process is known as the 
glass transition (GT), and although is very important from 
the fundamental and technological point of view, there are 
still many unsolved questions related to it Not all ma- 
terials are able to form glasses, and many criteria have been 
proposed in order to explain the ability of a material to reach 
the glassy state Q. The temperature where the GT occurs, 
is called the glass transition temperature (T g ). Many factors 
determine the T g , but among these, the chemical composition 
is fundamental. Chalcogenide glasses (formed with elements 
of the VI column) are very useful for understanding the ef- 
fects of the chemical composition J3j. In fact, T g can be raised 
or lowered by adding impurities, and the fragility of the glass 
can be changed from strong to fragile jij. Recently, by us- 
ing stochastic matrices the empirical modified Gibbs- 
DiMarzio law that accounts for the relation between T g and 
the concentration of modifiers jjj has been obtained, including 
the characteristic constant that appears in the law for almost 
any chalcogenide glass |1. 

The ease of glass formation in covalent glasses can be ex- 
plained at least in a qualitative way by the rigidity theory 
(RT), introduced first by Phillips bl and further refined by 
Thorpe ju^. By considering the covalent bonding as a me- 
chanical constraint, within this theory, the ease of glass for- 
mation is related with the proportion of available degrees of 
freedom and the number of constraints. If the number of con- 
straints is lower than the degrees of freedom, there are zero 



frequency vibrational modes called floppy |ll| , and the result- 
ing network is under-constrained. A transition occurs when 
the lattice becomes rigid, and at the corresponding chemical 
composition, the glass is easy to form. Many features of this 
transition have been experimentally observed Q Jl^]. Also, 
one of the authors proposed that rigidity can be related with 
the statistics of the phase space energy landscape jy^, since 
the number of floppy modes is equal to the number of differ- 
ent configurations of the system with nearly equal minimal 
energies fl3|, and thus is a way to evaluate the function that 
gives the number of minima energy basins |l4| . 

In a recent paper, Thorpe et. al. ju^], remarked that in real 
glasses, even though formed at relative high temperatures, 
where the entropic effects are dominant, it is not correct to 
fully ignore energetic contributions which can favor particu- 
lar structural arrangements over others, (e.g. in a binary sys- 
tem chemical aggregation between unlike particles favor local 
chemical aggregation) . One interesting question that they ad- 
dress is how the structure itself can incorporate non-random 
features in order to minimize the free energy at the temper- 
ature of formation. They answer this question by proposing 
that the structure can self-organize avoiding stress in the ran- 
dom formed network Jl5[ |. In the literature exists experimen- 
tal evidence for self-organization in glasses |H| , this evidence 
has been associated with the intermediate phase proposed by 
Phillips p7| . In a previous work pj] , we observed that in 
a model of an associative fluid (the Cummings-Stell model), 
some thermodynamics features can be associated with a rigid- 
ity transition, and in particular, it was shown that a glass 
transition occurs very near to the RT. Also, we showed that 
using the MC step as a time parameter in a NPT ensemble, 
we were able to control the cooling rate of a liquid melt in 
a qualitative way. In this work, we go further by looking at 
the self-organization properties in the Cummings-Stell model 
(CS), using a Monte-Carlo (MC) computer simulation in a 
grand canonical (GC) ensemble. Compared with the NPT 
ensemble, the GC ensemble has the advantage of reaching 
equilibrium faster |m| , given the opportunity to visit a wider 
range of equilibrium an non-equilibrium phases. In order to 
talk about the thermodynamic properties of these phases, we 
basis our assumptions in the fact that their life time is larger 
than the observation time (averaged time) This time is 

also larger than the molecular relaxation time, which we can 
adjust by tuning the MC steps of single particle movements, 
and the MC steps of the formed clusters movements. As a 
consequence, the slower a liquid is cooled, the longer the time 
available for configuration sampling at each temperature, and 
hence the access to the homogeneous nucleation which leads to 
crystallization. We point that this nucleation produces stress 
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in the obtained structure, as a counting of floppy modes re- 
veals. As a counter part, the faster the liquid is cooled, there 
is less time available for homogeneous nucleation and hence 
less stress is produced in the structure, inducing a local self- 
organization. This framework allow us to address the question 
of what are the structural and thermodynamic properties of a 
self-organized structures and how they behave as the system 
is cooled. With this in mind, we perform MC simulations 
where configurations that produce stress in the system are 
rejected, in a similar way to that proposed by Thorpe et. al. 
mM for studying rigidity self-organization. As we will see, our 
results are in agreement with ref. Jl5|, since the avoidance 
of stress, bias the system to a glass state. The layout of this 
work is as follows, in section II we introduce the model to be 
used, in section III a method for indirect controlling of the 
various relaxation times is introduced, and in section IV we 
discuss the effect of self-organization. Finally, in section V 
the conclusions of the work are given. 



II. MODEL 

We choose a simple model of an associative fluid: the 
Cummings-Stell model (CSM) of a two component system 
(A and B) of associating disks in 2D, both of the same size 
[|19|. The particles interact via a potential permitting core 
interpenetration of the A and B monomer discs, so that the 
bond length L is less than the core diameter a. Without loss 
of generality we assume a = 1. The interactions are given as 
follows: 



Uijir) = EC(r) + (1 - *«)U«.(r) 



U h A d A {r) 



U h / B {r) 



U as (r) 



U 



U h BA {r) = 




r < L — Q.5w, 
L - 0.5w < r < 1, 
r > 1, 



r < L - 0.5w, 

-E a3 - D L- 0.5w <r < L + 0.5w, 
r>L + 0.5w, 



where i and j stand for the species of the particles and take 
values A and B. r is the separation between centers, L is the 
bonding distance and w is the width of the attractive intra- 
core square well (fig. 1). The model allows the formation of 
dimer species for small values of the bonding length parame- 
ter, the formation of chains, if the bonding length is slightly 
larger, and also vulcanization with fixed maximum coordina- 
tion number for different bonding length values close to the 
diameter of particles, as shown in fig. 1. In order to be able 
to fix a maximum coordination number in each simulation , 
we take D — * oo as was done before in other works ]l9fl , [ p2| , 
[ p3| . This choice has the effect that unlike particles avoid 
bond-lengths between L + 0.5w and 1, and thus coordinations 
higher than a desired maximum are not allowed. Numerically, 
this condition means that in the MC simulations, we never 
consider bond distances in the previous range. The corre- 
sponding values for each maximum coordination number are 
given in the following table, 



^max 


L 


w 


cpx 3 


0.65 


0.1 


cpx 4 


0.78 


0.1 


cpx 5 


0.91 


0.1 



Table 1. Parameters of the CS model that fix the maxi- 
mum coordination of the particles as used in this work. The 
notation cpx r, means complex of particles with maximum 
coordination r max . 



III. GLASS TRANSITION BY CONTROLLING 
RELAXATION TIMES 

We start by pointing out that a supercooled liquid phase 
is metastable with respect to the crystalline state, and this 
supercooling can be achieved if nucleation is inhibited during 
cooling |HJ . One way of inhibit nucleation is by perform- 
ing a rapid quench of the liquid; in these terms, two char- 
acteristic times n (the time required for crystallization) and 
T2 (the time corresponding to molecular relaxation) compete 
between crystallization and vitrification j2lj. In this work, we 
simulate this effect in two different ways, which we will show 
that at the end turn out to be very similar: one is to use 
the MC steps as a time parameter in the GC ensemble [ pij , 
where we tune the ratio between n and T2 in an indirect way, 
by controlling the ratio between steps of particle and cluster 
rearrangements, since the first is the most important factor 
for molecular relaxation, while the second optimize crystal- 
lization. The second way which we observed that leads to 
supercooling is the self-organization of rigidity, as we will see 
in the next section. 

To implement supercooling using the MC steps by indi- 
rect control of the relaxation times, we use a GC Metropo- 
lis Monte-Carlo method ^5|. The procedure has two nested 
loops. In the inner one, the particles are moved inside the vol- 
ume, and an interchange of particles with the particle reser- 
voir is allowed. This loop is performed N m times. The particle 
movements inside the volume allow to rearrange the struc- 
ture, and thus this is related with the molecular relaxation 
of the structure (ra). In the outer loop, cluster rearrange- 
ments and the average of the thermodynamical quantities are 
performed, each time that N m cycles of the inner loop are fin- 
ished. The external loop is related with crystallization, since 
cluster movements promote the growth of bigger clusters. Is 
clear that if N m is high enough, the probability of having 
local cluster nucleation is high, and thus cluster movements 
allow to form a crystal by successive aggregation of small clus- 
ters. When N m is small, the local configurations are not in 
equilibrium, and hence the cluster movements promote the 
generation of a random network. 

In fig. 2, we show the results of the inverse density (p _1 ) 
against reduced temperature (T* = (/3e as ) _1 ) for several N m 
cycles which simulate different T2 times. In this figure, we 
fixed the parameters of the CSM that allows as maximum 
coordination four (cpx4), L = 0.78, w = 0.1, restricted to 
equimolar concentration !3\ia = fiflB ~ —0.4. Finally, af- 
ter each N m steps of the inner cycle, we allow the possibility 
of one cycle cluster rearrangement. In order to simulate the 
same cooling rate with several T2 times, we fix the outer cy- 
cle for 100 steps during the quasi-equilibration run, and 2000 
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times for a productive run. In such a way, we averaged over 
the same number of configurations for each different N m . As 
can be seen, for N m — 1500 a crystallization is observed, while 
for the other values, a glass transition is obtained, as is re- 
vealed by the figure and by a direct inspection of the resulting 
structures. In all the procedure, MC steps are controlled to 
have an acceptance ratio between 20 — 30%. 

An interesting observation, is that the fraction of particles 
with maximum coordination X4, depends strongly on N m . In 
table 2, we show this fraction against temperature for the 
same simulations presented in fig. 2. When the crystal is 
formed, \i is nearly one, while for the supercooled liquid at 
the same temperature, \i remains at the same order of magni- 
tude. We can understand this effect as follows: if full thermal 
equilibrium of the system is not allowed, is not possible to 
access the global minimum of the energy potential (2^], and 
hence the nucleation is prohibited, with the consequence that 
the supercooled liquid is structurally arrested at a finite tem- 
perature and restricted to explore the configurational space 
correspondent to a single basin in the energy landscape. In 
the next section, we use this idea to relate this observation 
with the self-organization of rigidity. 



T* 


N m = 12 


N m =40 


N m 


= 150 


N m = 1500 


0.40 





0011 


0.00113 





0025 


0.0019 


0.38 





0046 


0.00103 





0053 


0.0040 


0.36 





0022 


0.00404 





0043 


0.0079 


0.34 





0007 


0.01287 





0152 


0.0188 


0.32 





0077 


0.01677 





0495 


0.0591 


0.30 





0108 


0.04376 





1353 


0.9395 


0.28 





0258 


0.05391 





1538 


0.9786 


0.26 





0366 


0.05138 





2229 


0.9666 


0.24 





0725 


0.07417 





3100 


0.9487 


0.22 





1057 


0.08355 





3292 


0.9622 


0.20 





1116 


0.07762 





3496 


0.9617 



Table 2. Fraction of maximum coordinated particles as 
a function of the reduced temperature (T*) and molecular 
relaxation time, controlled by the parameter Nm. 



counting, goes as follows: since each of the r bonds in a site of 
coordination r is shared by two sites, there are r/2 constraints 
due to distance fixing between neighbours. If the angles are 
also rigid, in 2D there are (r — 1) constraints, to give, 



/ 



2iV - Nc 
2N 



< r > 



(r - V)x r 



where the last term corresponds to the angular constraints, 
Xr is the fraction of particles with coordination r, and < r > 
is the average coordination number, defined as, 



< r > = 



A rigidity transition occurs when / = and the system pass 
from a floppy network to rigid one. If / is a negative number, 
i.e., if there are more constraints than degrees of freedom, 
the lattice is overconstrained and / is the number of stressed 
bonds. In 2D, the rigidity transition leads to the critical 
value < r >= 2.0 if all angular constraints are considered, 
and < r >= 4.0 if the angular restoring forces are not strong. 

Within the Cummings-Stell model, rigidity comes from the 
association of particles: each bond generates a constraint, and 
the angular constraints are only produced by geometrical hin- 
drance, i.e., the angles between particles can change without 
a cost in energy, but within certain limits imposed by the re- 
striction of the hard-core interaction between like particles, as 
shown in fig. 3. For maximum coordination four, this means 
that only sites with coordination four have a contribution to 
angular constraints. It is true that sites with coordination 
two and three in principle should provide extra angular con- 
straints, since the hard core interaction gives a minimum angle 
between particles. However, the angles are not fixed and they 
have a wide region to allow particle movements, and thus do 
not contribute to the restriction counting. 

Taking into account the geometrical hindrance of the 
model, the number of floppy modes is now given by, 



/ 



< r > 



E 



<5rr max (r - 1)3 



IV. GLASS TRANSITION BY 
SELF-ORGANIZATION OF RIGIDITY 

The fact that the maximum coordination is not achieved 
for most of the particles in the supercooled liquid, means that 
inhibit nucleation is a natural way of inhibit crystallization, 
as was discussed in the introduction. This simple idea can 
be put in contact with the rigidity ideas of Phillips |nj and 
Thorpe jyj. As we mention before, in this theory, the ability 
for making a glass is optimized when the number of freedom 
degrees, in this case 2N, where N is the number of particles, 
is equal to the number of mechanical constraints (N c ) that 
are given by the bond length and angles between bonds. 

(2N — N C )/2N gives the fraction of cyclic variables of the 
Hamiltonian, and also corresponds to the number of vibra- 
tional modes with zero frequency (/), called floppy modes, 
with respect to the total number of vibrational modes. The 
counting of floppy modes in a mean- field, known as Maxwell 



where r max is the maximum allowed coordination, and 
<5 rT - max is a Kronecker delta. From here, is clear that when 
in a cluster we have a site with maximum coordination, rigid- 
ity raises since the delta function is different from zero and 
more constraints are added that over-constraint the cluster. 
Thus, sites with maximum coordination nucleate rigidity and 
produce stress in the lattice. For example, in a crystal with 
coordination four, ^4 = 1 an d / = —3/2, which means that 
the lattice is over-constrained. According to Phillips, when 
/ = 0, it is easy to form a glass, since the material is neither 
over-constrained (that produce explosive exothermic crystal- 
lization due to strain energy |^|) nor under-constrained (lead- 
ing to the formation of a molecular crystal), the system is 
trapped in a configurational limbo where fluctuations do 
not provide a pathway to the crystalline phase. 

Furthermore, in the last section we have showed that the 
fraction of particles with maximum coordination is in close 
connection with the molecular relaxation time, which in other 
words means that to form a glass, nucleation of stress must 
be prevented. From our previous results, we can observe that 
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the probability of formation of a nucleated structure is small 
due to the high molecular relaxation time. We decided to 
follow these ideas by proceeding in the opposite way than in 
the last section, i.e., we inhibit rigidity nucleation by rejecting 
configurations with maximum coordination and then we see if 
we are able to bias the system to a glassy state, in such a way 
that we simulate long molecular relaxation times. Observe 
that rejecting configurations that produce stress is the same 
process of self-organization that was considered by Thorpe et. 
al. in order to form stress-free lattices jlB) . In that sense, we 
look if self-organization of rigidity is able to produce a glass 
transition. This kind of simulation is usually called biased 
Monte-Carlo Q. 

To study the effect of self-organization, we made the same 
MC procedure described in the previous section, but with 
N m fixed to the value that gives crystallization (N m — 1500). 
The only difference with the previous case is that now we 
reject particle movements that produce a site with maximum 
coordination. 

In figure 4, we present the behaviour of the inverse of 
the density as a function of the reduced temperature 

T* = (Peas)" 1 , with the condition that allows maximum 
four neighbours (cpx4), restricted to equimolar concentration 
fi[iA = PfJ-B = —0.4, (open squares). As the temperature 
is slowed down, we can observe a continuous decrease in p~ . 
However, for reduced temperatures lower that 0.30, a jump in 
p _1 is observed when all the configurations are allowed. This 
jump corresponds to the crystalline like phase transition, as 
can be argued by the shape of the transition, from an inspec- 
tion of the configuration obtained, and by the radial distri- 
bution function. Due to the fact that it is possible to keep 
the system without stress, we develop the same simulation 
as before but rejecting in the simulation every configuration 
that contains a particle with coordination four. The results 
are presented in fig. 4 with dashed squares. In that case, the 
system remains as a supercooled fluid. Moreover, the system 
can not form a crystal structure as occurs in the simulation 
in which we allow stress, and thus do not present a usual 
phase transition, instead a glass like transition is observed. 
These results shows that self-organization of rigidity is able 
to produce a glass. 

An important remark is that avoiding configurations with 
maximum coordination is not equivalent to consider a CSM 
without self-organization but with a lower maximum coordi- 
nation. For example, in fig. 4 we plot the results of a sim- 
ulation without rejection for a CSM that allows maximum 
coordination three (cpx3). As can be seen, the model also 
presents crystallization. 

In figure 5, we present the results for the same kind of 
simulation but for a system that allows maximum coordina- 
tion three (cpx3). As can be seen, the rejection of stressed 
configurations also leads to a glass like transition. 

Now we turn our interest to the condition of maximum co- 
ordination five (cpx5) , as show in fig. 6. As can be seen, in this 
case the glass transition is produced even when the stressed 
configurations are rejected. This fact can be understood in 
terms of rigidity in the following way: when r max = 5, it 
is impossible to have a crystal due to geometric frustration 
at equimolar condition = P^b = —0.4, and xs "C 1. 

Since the rigidity transition without angular restrictions oc- 



curs only when < r >= 4, most of the configurations do not 
produce stress and the system behaves freely (i.e. we do not 
need to reject any configuration) as a self-organized system. 
In this sense, geometric frustration induce self-organization of 
the system. As a corroboration of this fact, in fig. 7 we show 
the number of floppy modes as a function of the average co- 
ordination number, using the Maxwell counting. We remark 
that each coordination corresponds to a certain temperature 
of the simulation. For example, at high temperatures, all the 
models with different maximum coordination fall in the same 
fine, since in the liquid the probability of nucleation is very 
low. However, for the case of cpx5, cpx3 and cpx4 without 
stress, all the simulations fall again in the same line even for 
low temperatures, since the self-organization means that the 
clusters grow without angular constraints (stress free) . When 
this line is extrapolated to / = 0, we obtain < r >= 4, 
which is the value for a rigidity transition without angular 
constraints. If angular constraints are allowed, the simula- 
tions for low temperatures falls outside the line determined by 
self-organization, and the rigidity transition occurs at lower 
values of < r > . Finally, we can compare these results to the 
floppy mode counting made for the glass transition using the 
method of tuning the different relaxation times. In fig. 8, we 
present the number of floppy modes as a function of < r >. 
As can be seen, when N m is high, there is a transition of 
rigidity due to nucleation, while for low N m the system tends 
to stay in the line of self-organization. 

V. CONCLUSIONS 



In this work, we have explored the connection between self- 
organization of rigidity, and the supercooling of a liquid to 
form a glass. By considering an associative fluid model, we 
showed that the competition between two different character- 
istic times, molecular relaxation and crystallization times, can 
be modelled using a MC simulation, where the number of cy- 
cles between particle and cluster moves is controlled. The re- 
sults of these simulations, suggested that nucleation is avoided 
during supercooling and produce self-organization in the sense 
of the rigidity theory. This idea was also tested by making 
MC simulations but avoiding stressed configurations. As a 
result, we were able to produce glass transitions using self- 
organization. In a model (cpx5) with geometrical frustration, 
this self-organization is provided by geometry, and thus glass 
transition occurs without rejecting configurations. All of the 
results of this article are in agreement with the Phillip 's idea 
that glass transition is related with rigidity due to the lack of 
a pathway to crystallization j27|. Many of these facts, can 
also be studied form an energy landscape point of view, as we 
will show in future works. 
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l.-Cummings- Stell model. 



2. - Inverse of the density (p -1 ) as a function of the scaled 
temperature (T*) for different values of N m . 

3. -Counting of angular constraints in the Cummings-Stell 
model. A cluster of two particles has no angular constraints, 
since one of the particles can rotate 360° around the other, 
while a in a cluster with coordination four the angle between 
particles is fixed, which leads to 3 angular constraints. 

4. -The same as fig. 2 for the case of stress free nucleation 
for a system with maximum allowed coordination four (cpx4) , 
with and without stress (squares and dashed squares). We 
include a simulation for maximum coordination three (cpx3) 
with stress. 

5. - Inverse of the density (p -1 ) as a function of the scaled 
temperature (T*) for maximum coordination three with and 
without stress (triangles and dashed triangles). 

6. -Inverse of the density (p^ 1 ) as a function of the scaled 
temperature (T*) for maximum coordination five (cpx5) with 
and without stress (pentagons and dashed pentagons). For 
comparison purposes, we include a simulation of cpx4 without 
rejecting configurations (squares). 

7. -Number of floppy modes as a function of the coordina- 
tion number (< r >) for each of the models with different 
maximum coordination number. 

8. -Number of floppy modes as a function of < r > for sev- 
eral N m relaxation times without rejecting any configuration. 
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